home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
L' Effet Pommier 3
/
L'Effet Pommier - Volume 03.iso
/
Programmation
/
gray image 2.1
/
median_filter.cc
< prev
next >
Wrap
Text File
|
1995-05-31
|
9KB
|
297 lines
// This may look like C code, but it is really -*- C++ -*-
/*
************************************************************************
*
* Grayscale Image
*
* Perform median filtration of an image
*
* Implementations of several flavors of median filtration: by rows
* only, filter only columns, filter rows and then filter columns. The
* latter is equivalent to a 2D median filtration with a diamond-shaped
* window. 2D filtration with a square window fits the existing framework,
* too, but not yet implemented.
*
* When the filtration window sticks out of the image, we extrapolate
* the corresponding pixels. Incidentally, it means that the border
* pixel(s) are never affected by the filtration. Indeed, say, for a
* 3-point row filtration
* new_pixel[i] = median_of(pixel[i-1],pixel[i],pixel[i+1])
* if we assume that pixel[-1] = pixel[0], then
* median_of(pixel[0],pixel[0],pixel[1]) is always pixel[0]
* (majority wins). The same holds for other types of filtration.
*
* Square window 2D filtration:
* For every pixel of the source image x[i,j] compute the
* pixel y[i,j] of the filtrated image according to the
* rule
* y[i,j] = median_of( x[i+p,j+l] ), p,l = -ws...0...ws
*
* where 2*ws+1 is the window size. The median of the sequence
* is the middle term of the ordered sequence (in our case the
* sequence always contains the odd number of terms).
* This isn't implemented yet.
*
* $Id: median_filter.cc,v 2.5 1995/03/14 22:11:03 oleg Exp oleg $
*
************************************************************************
*/
#include "filter.h"
#include "minmax.h"
/*
*------------------------------------------------------------------------
* Simple 3-point median filtration
*/
// Find a median of three numbers
static inline GRAY_SIGNED median_of(const GRAY_SIGNED a1,
const GRAY_SIGNED a2,
const GRAY_SIGNED a3)
{
if( a3 >= a1 )
if( a2 >= a3 )
return a3; // order is a1, a3, a2
else
if( a2 > a1 )
return a2; // order is a1, a2, a3
else
return a1; // order is a2, a1, a3
else
if( a2 >= a1 )
return a1; // order is a3, a1, a2
else
if( a2 > a3 )
return a2; // order is a3, a2, a1
else
return a3; // order is a2, a3, a1
}
// 3-point row-wise filtration
IMAGE& FilterIt::median_3_row(void)
{
if( image.ncols < 2 )
image.info(),
_error("we need at least 3 columns for 3-point median row filtration");
register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
register GRAY_SIGNED prev_pixel; // in the row
for(GRAY_SIGNED * row_end = (GRAY_SIGNED *)image.pixels + image.ncols;
row_end <= (GRAY_SIGNED *)image.pixels + image.npixels;
row_end += image.ncols)
{
prev_pixel = *pp++; // At the beginning, pixel[-1]=pixel[0]
// and the first pixel never changes
while(pp < row_end-1)
{
GRAY_SIGNED new_pixel = median_of(prev_pixel,pp[0],pp[1]);
prev_pixel = pp[0];
*pp++ = new_pixel;
}
pp++; // if we assume p[ncols] = p[ncols-1]
// then the last pixel is always the
// same
}
assert( pp == (GRAY_SIGNED *)image.pixels + image.npixels );
return image;
}
// 3-point col-wise filtration
IMAGE& FilterIt::median_3_col(void)
{
if( image.nrows < 2 )
image.info(),
_error("we need at least 3 rows for 3-point median col filtration");
register const int ncols = image.ncols; // caching
const int fl_jump = image.npixels - ncols; // between 1. and last rows
register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
register GRAY_SIGNED prev_pixel; // in the col
GRAY_SIGNED * last_row = (GRAY_SIGNED *)image.pixels + fl_jump;
for(; last_row < (GRAY_SIGNED *)image.pixels + image.npixels; last_row++)
{
prev_pixel = *pp, pp += ncols; // At the beginning, pixel[-1]=pixel[0]
// and the first pixel always stays
for(; pp < last_row; pp += ncols)
{
GRAY_SIGNED new_pixel = median_of(prev_pixel,pp[0],pp[ncols]);
prev_pixel = pp[0];
pp[0] = new_pixel;
}
pp -= fl_jump-1; // if we assume p[nrows] = p[nrows-1]
// then the last pixel is always the
// same
}
assert( pp == last_row - fl_jump );
return image;
}
/*
*------------------------------------------------------------------------
* Simple 5-point median filtration
*
* Note about the algorithm:
* we can always reduce finding a median of 5 numbers to finding a median
* of three numbers by noting that max and min of four numbers never
* count. Indeed, if mx = max(a1,a2,a4,a5) and a3 > mx, then
* the order is "? ? ? mx a3" and mx won't figure in in the median, so
* the median is the same as of "? ? ? a3". If a3 < mx, then the order
* is "? ? ? a3 mx" and mx won't count even more. The same holds for the
* min of four numbers.
*/
// Find a median of 5 numbers
static inline GRAY_SIGNED median_of(const GRAY_SIGNED a1,
const GRAY_SIGNED a2,
const GRAY_SIGNED a3,
const GRAY_SIGNED a4,
const GRAY_SIGNED a5)
{
register GRAY_SIGNED mn1 = a1, mx1 = a2;
if( a1 > a2 )
mn1 = a2, mx1 = a1;
register GRAY_SIGNED mn2 = a4, mx2 = a5;
if( a4 > a5 )
mn2 = a5, mx2 = a4;
return median_of(min(mx1,mx2),a3,max(mn1,mn2));
}
// 5-point row-wise filtration
IMAGE& FilterIt::median_5_row(void)
{
if( image.ncols < 4 )
image.info(),
_error("we need at least 5 columns for 5-point median row filtration");
register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
register GRAY_SIGNED prev_pixel1; // pixel[*,j-1]
register GRAY_SIGNED prev_pixel2; // pixel[*,j-2]
for(GRAY_SIGNED * row_end = (GRAY_SIGNED *)image.pixels + image.ncols;
row_end <= (GRAY_SIGNED *)image.pixels + image.npixels;
row_end += image.ncols)
{
prev_pixel2 = prev_pixel1 = *pp++; // pixel[-2] = pixel[-1] = pixel[0]
// and the first pixel always stays
while(pp < row_end-2)
{
GRAY_SIGNED new_pixel =
median_of(prev_pixel2,prev_pixel1,pp[0],pp[1],pp[2]);
prev_pixel2 = prev_pixel1;
prev_pixel1 = pp[0];
*pp++ = new_pixel;
}
*pp = median_of(prev_pixel2,prev_pixel1,pp[0],pp[1],pp[1]);
pp += 2; // pp[ncols-1] is unchanged
}
assert( pp == (GRAY_SIGNED *)image.pixels + image.npixels );
return image;
}
// 5-point col-wise filtration
IMAGE& FilterIt::median_5_col(void)
{
if( image.nrows < 4 )
image.info(),
_error("we need at least 5 rows for 5-point median col filtration");
register const int ncols = image.ncols; // caching
const int fl_jump = image.npixels - 2*ncols; // between 1. and last b1 rows
register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
register GRAY_SIGNED prev_pixel1; // pixel[i-1,*]
register GRAY_SIGNED prev_pixel2; // pixel[i-2,*]
register GRAY_SIGNED next_pixel1; // pixel[i+1,*]
register GRAY_SIGNED next_pixel2; // pixel[i+2,*]
GRAY_SIGNED * last_b1_row = (GRAY_SIGNED *)image.pixels + fl_jump;
for(; last_b1_row < (GRAY_SIGNED *)image.pixels + image.npixels - ncols;
last_b1_row++)
{
prev_pixel1 = prev_pixel2 = *pp, pp += ncols; // skip over pixel[0,*]
// which is unaffected
next_pixel1 = pp[ncols];
for(; pp < last_b1_row; pp += ncols)
{
next_pixel2 = pp[2*ncols];
GRAY_SIGNED new_pixel =
median_of(prev_pixel2,prev_pixel1,pp[0],next_pixel1,next_pixel2);
prev_pixel2 = prev_pixel1;
prev_pixel1 = pp[0];
next_pixel1 = next_pixel2;
pp[0] = new_pixel;
}
pp[0] = median_of(prev_pixel2,prev_pixel1,pp[0],next_pixel1,next_pixel1);
pp -= fl_jump-1; // skip over pixel[nrows-1]
}
assert( pp == last_b1_row - fl_jump );
return image;
}
/*
*------------------------------------------------------------------------
* Median Filtration drivers
* pick up and execute the most efficient method given input parameters
*/
IMAGE& FilterIt::median(const Direction how, const int win_size)
{
image.is_valid();
switch(win_size)
{
case 3:
switch(how)
{
case RowsAndColumns:
return median_3_row(), median_3_col();
case Columns:
return median_3_col();
case Rows:
return median_3_row();
};
break;
case 5:
switch(how)
{
case RowsAndColumns:
return median_5_row(), median_5_col();
case Columns:
return median_5_col();
case Rows:
return median_5_row();
};
break;
default:
_error("Sorry, not implemented yet. "
"FYI: implemented windows: 3 and 5, and they are fast");
}
_error("There is no way we can reach here...");
return image; // just to make the compiler happy
}
#if 0
/*
*------------------------------------------------------------------------
* Finding a median of the sequence of GRAYs
*
* Algorithm
* Described in the book "Numerical Recepies".
* The algorithm does not require sorting the sequence. The main idea is
* that if the median of the sequence is a "middle" term of the sequence,
* the following relation should hold
*
* SUM{ (xmed - x[i]) / |xmed - x[i]| } = 0
*
*/
static GRAY find_median(GRAY seq[], const int no_elems)
{
}
#endif